Short-distance properties of Coulomb systems 
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We use the operator product expansion to derive exact results for the momentum distribution 
and the static structure factor at high momentum for a jellium model of electrons in both two and 
three dimensions. It is shown that independent of the precise state of the Coulomb system and for 
arbitrary temperatures, the asymptotic behavior is a power law in the momentum, whose strength 
is determined by the contact value of the pair distribution function g{0). The power-law tails are 
quantum effects which vanish in the classical limit ft — > 0. A leading order virial expansion shows 
that the classical and the high-temperature limit do not commute. 
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I. INTRODUCTION 



The basic constituents of ordinary matter are electrons 
and nuclei combined in such a way that there is no net 
overall charge. Within a non-relativistic approximation 
and treating the nuclei as point particles, the interaction 
is fully described by an instantaneous Coulomb potential 
^ e^/r at arbitrary distances. In spite of the long-range 
nature of this interaction and the divergent attractive 
force between electrons and nuclei at short distances, one 
expects an overall neutral Coulomb system to be stable 
in the sense that the ground state energy (or free energy 
at finite temperature) is finite and scales linearly with the 
total number N of particles. It is one of the major accom- 
plishments of theoretical physics to show that - beyond 
the exactly solvable case of the hydrogen atom - these 
expectations can indeed be proven rigorously. The proof 
crucially relies on the fact that electrons are fermions and 
are thus constrained by the Pauli principle^'. Since nei- 
ther the size and mass nor the statistics of the nuclei play 
a role in this context, a simple approximation which cap- 
tures much of the basic physics of Coulomb systems is the 
well-known jellium model, where the nuclei are treated as 
a homogeneous background that precisely cancel the neg- 
ative charge of the Coulomb gas of electrona^l. At zero 
temperature, this model is fully specified by the stan- 
dard dimensionless interaction strength r^ = ro/ao- Here 
ao = f? jme^ is the Bohr radius while rg is the average 
spacing between electrons, connected with the electron 
density n via tq = (3/47rn)^/'^ in 3D and rg = (l/7m)^/^ 
in 2D, respectively. Despite the fact that the periodic 
arrangement of the nuclei is ignored, the jellium model 
provides a reasonable starting point to describe elemen- 
tary properties of metals like their cohesion energy or 
compressibilitySSl^ Unfortunately, however, beyond the 
fundamental issue of stability and extensivity, there are 
hardly any exact results even for this highly simplified 
model. It is only in the high-density limit r^ ^ 1 where a 
perturbative expansion around the non-interacting Fermi 
gas is possible. A simple argument for this is provided 



by writing the jellium Hamiltonian 



H = -- 



Hy, 



i<2 



(1) 



in dimensionless form, with Ry= e^ /2aQ as the unit of en- 
ergy and particle coordinates r^ measured in units of tq 
(both the energy iJb of the background as well as the in- 
teraction energy between electrons and the background 
are constants and thus need not be written explicitly). 
Clearly, as r^ — >■ 0, the kinetic energy dominates and the 
Coulomb interaction can be treated within perturbation 
theory. The expected ground state is a Fermi liquid, with 
a finite jump < Z < 1 of the momentum distribution at 
a spherical Fermi surface |k| = kp^. While typical values 
of Ts w 1 — 5 in metals^l are outside the range of perturba- 
tion theory, at least the qualitative features of electrons 
in metals are captured correctly in this picture. For very 
large values of r^, the uniform electron liquid is expected 
to eventually form a Wigner crystal, which minimizes the 
interaction energy in (fl]). In addition, non-trivial phases 
such as anisotropic quantum liquid crystals are likely to 
appear at intermediate values of r^. Indeed, in two di- 
mensions a direct transition from a uniform electron liq- 
uid to a Wigner crystal as a function of r^ can be ruled 
out by a quite general thermodynamic argumenlP. More- 
over, even in the Wigner crystal, the electron spin gives 
rise to strong quantum fluctuations due to ring exchange 
processes, leading to a complex magnetic structure, see 
e.g. Refs.EHm. 

In view of the still poorly understood phase diagram 
of even the simplest model for a many-body system with 
Coulomb interactions, it is of considerable interest to 
derive exact relations that hold independent of the in- 
teraction strength and the particular state in question. 
Our aim in the present work is to show that such re- 
lations follow rigorously from the operator product ex- 
pansion (OPE) of quantum field theory. They constrain 
the short-distance physics of Coulomb systems in a man- 
ner which is analogous to the so-called Tan relational^Hlll 
for fermions with interactions that have effectively zero 



range. Physically, these relations rely on the fact that 
at sufficiently short distances only the two-body interac- 
tion is relevant. At the level of the many-body ground 
state wave function for Coulomb systems, this observa- 
tion goes back to KimbalP^ti^l, Beyond providing a clear 
mathematical proof of these relations, a novel and quite 
non-trivial aspect of our present derivation via the OPE 
is the fact that the relations will be shown to be valid 
at the level of operators. As a result, they apply to any 
state of the system, e.g., to a Fermi liquid or a Wigner 
crystal, at zero or at finite temperature and also in a 
few-body situation. The only change is the value of the 
'contact' g(0), a dimensionless measure of the probabil- 
ity that two electrons with opposite spin are found at a 
coincident point in space. 

In order to illustrate the universal features in the short- 
distance behavior of Coulomb systems, it is instructive 
to consider the hydrogen atom as a simple and exactly 
solvable system. In fact, the basic relations that will 
subsequently be proven for the many-body case show up 
already in this elementary textbook problenJi^. The wave 
function for relative motion in the hydrogen atom has the 
well-known forirPSI 



V'nim(r) 
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where Y/ml^, <p) and L^_j^_^{x) are spherical harmonics 
and generalized Laguerre polynomials, respectively (we 
assume the proton to be infinitely heavy, so the reduced 
mass is equal to the electron mass m). Its Fourier trans- 
form Tnim (q) has been calculated by Podolsky and Paul- 
ing2Sl. It is given by 

T„,™(q) ^2^'+V(aon)^/^(~»)'/!j ^^^"~^^,^^ ^l^„.(^, <^) 



C ^i+i 



{e + if+^ ""'-'Vc^ + v' 



C^-l 



(3) 



where C, = qaon, and C^i_-^{x) denotes a Gegenbauer 
polynomial. The momentum distribution of the elec- 
tron is the absolute square of the momentum space wave 
function. Using Eq. ([s]), |T„/„i(q)p turns out to de- 
crease asymptotically as l/g*+^' for large momentum 
C = qaon^ 1. The leading order term 



|T„oo(q) 



IV-nOoWl^ 



+ Oil/q'') (4) 



in the momentum distribution therefore only involves the 
contribution from s-states. They are the only ones with a 
finite probability density |'!/'nOo(0)P — l/7r(aon)^ for the 
electron and proton to be found at a coincident point in 
space. 



Remarkably, the same 'contact' density also appears 
in the high-momentum tail of the atomic form factor 
Pnim (q) , which is the Fourier transform of the electronic 
density distribution. Its leading contribution at large q 



Pnoolqj = I d-^re ''»■'" IV'nOo(r) I 
16tt |V'nOo(0)P 



ao 



+ 0{l/q') 



(5) 



comes again from s-states, while higher angular momenta 
are associated with faster decaying power laws. As for 
the momentum distribution, the coefficient of the high- 
momentum tail contains the contact density |V'nOo(0)P- 
Moreover, both the momentum distribution and the form 
factor depend only on the magnitude q = |q| of the wave 
vector, i.e., they have spherically symmetric tails since 
only s-states contribute. As will be shown in the follow- 
ing, the power laws found in the hydrogen atom and the 
fact that the physics at short distances is rotation in- 
variant also show up in the many-body context, even for 
inhomogeneous or anisotropic phases. More precisely, the 
momentum distribution is replaced by the Fourier trans- 
form of the one-particle density matrix, while the atomic 
form factor becomes the static structure factor S'(q) of 
the many-body system. 

The article is structured as follows: in Sec. |Tlj we in- 
troduce the jellium model, the one- and two-particle den- 
sity matrix as well as some basics of the operator prod- 
uct expansion. In Sec. |III[ the OPE is used to derive 
the exact short-distance behavior of general many-body 
wave functions. Moreover, it is shown that this implies 
power-law tails in both the momentum distribution and 
the static structure factor which depend on the particular 
state in question only through the value of the contact 
g{0). A direct computation of the short-distance OPE of 
the density-density correlator and the one-particle den- 
sity matrix is presented in Sec. |IV[ Finally, an explicit 
calculation of the contact in both the classical and the 
high-temperature limit is given in Sec.lV] It is shown that 
for a Coulomb system these limits give opposite results 
and thus do not commute. The article is concluded by a 
summary and outlook, Sec. |VI| The Feynman rules of the 
diagrammatic calculation and some details of the evalu- 
ation of some few-particle matrix elements are discussed 
in Apps.[^and|B] 



II. JELLIUM AND OPE 

In second-quantized form, the Hamiltonian of the jel- 
lium model is given by 

2y72 



H = H^ + j d'^x ^t _|^^^(x) 



(6) 



where a suniination over spin indices a ^ti 4- is implied. 
Since the Hamiltonian involves only one- and two-body 
interactions, the expectation value of the energy in a state 
described by an arbitrary A^-body density matrix only in- 
volves the reduced one- and two-particle density matrices 
7^"'^^ and 7^^''^. In a spin-resolved form and in a position 
space representation, the former can be expressed as 



7^^^ I 



X, X 



(^l(x)V'.(x')> 



(7) 



Its Fourier transform with respect to x— x' then gives the 
momentum distribution (see Eq. ( |30p below). Regarding 
the two-particle density matrix 7'^!, 1'; 2, 2'), one only 
needs the diagonal elements 1 = 1', 2 = 2', which define a 
dimensionless, spin-resolved pair distribution function 

ri,(x)n„,(x').g^,.,(x,x') = (^t (x)^t^(x')^^,(x')V,^(x)) . 

(8) 
The pair distribution function is a measure of the prob- 
ability to find an electron with spin projection a' at 
position x' given an electron with spin projection a is 
at x. For a completely uncorrelated system one has 
(/(x, x') = 1. Note that there is no assumption here 
about translation invariance, which is certainly broken 
in the Wigner crystal. Taking into account the possibil- 
ity of a non-vanishing spin polarization C, = {n^ — n^)/n, 
the total pair distribution function 

5(x,x') =1 (1 + 0' gn(x,x') + i (1 - 0' .9u(x,x') 



+ ^(i-C')m(x,x') 



(9) 



of a spin one half Fermi gas is a weighted sum of contri- 
butions g^ ^, . They all approach unity as x — x' — )• cx3 
and so does g(x, x'). For small separations x — x' — >■ 0, in 
turn, the pair distribution function for equal spins van- 
ishes quadratically because of the Pauli principle. As a 
result, one finds 



5(0) 



^C 



C) m(o) 



(10) 



for the total pair distribution function at vanishing sepa- 
ration. Note that in a situation where the electronic state 
is not translation invariant, the local value g(Q) = ^(R, 0) 
of the pair distribution function depends also on the 
'center-of-mass' coordinate R = (x-|-x')/2, a dependence 
which is suppressed in the following. 

Both the one-particle density matrix and the pair dis- 
tribution function can be expressed as expectation values 
of operators at different points in space. The operator 
product expansion - specified here to the relevant case of 
equal times - provides an expansion of an operator prod- 
uct Oa Ob at nearby points in space in terms of local 
operators: 

a(R - |)a(R + |) = E W^«(r)0„(R). (11) 



It is important to emphasize that Eq. (Ill is an op- 
erator relation, i.e., it is valid for expectation values 
between any state. The state-independent coefficients 
M^„(r) are ordinary c- numbers and are called the Wilson 
coefficients. They depend both on n and the specific op- 
erators Oa and Ob which appear on the left-hand side 



of Eq. (11). The scaling dimension A„ of a local oper- 
ator On is defined by the property that the correlation 
between ©„ and its hermitian conjugate at points sep- 
arated by a small distance r and time t asymptotically 
behaves like t~^". For example, the operator ip^^ipcr has 
scaling dimension A = d. The values of A„ determine 
the dependence of the Wilson coefficients at small sepa- 
ration r via 



W„(r) 



„A„-A„-A6 



fir/ao) 



(12) 



where / is a function of the dimensionless ratio r/ao. The 
operators On with the lowest scaling dimension therefore 
govern the behavior of an operator product at small sep- 
aration. In particular, Wilson coefficients which are non- 
analytic in r give rise to power law tails of the associated 
correlator OaOb at large momentum. 

In practice, the Wilson coefficients may be determined 
by performing few-particle calculations. Indeed, since 
they are state-independent, it is sufficient to calculate the 
matrix element of Eq. ( 11 ) between simple (few-particle) 
states for which (On) 7^ 0. The coefficients W„(r) then 
follow by matching both sides of Eq. (11). As will be 
shown below, an operator of particular interest in the 
present context is the two-particle operator 



a(R) = 44'14'ii'ti'^) ■ 



(13) 



In analogy to the notion used for fermions with short 
range interactions, we shall refer to this as the contact 
operator. In contrast to the case where the interactions 
have zero range and Oc acquires an anomalous dimension 
twcpil, it has a finite expectation value in the presence of 
Coulomb interactions. As a result, the scaling dimen- 
sion Ac = 2d of the contact operator is the one inferred 
from simple dimensional analysis. The contact is a cen- 
tral quantity which determines the leading short-distance 
singularities of Coulomb systems and, in particular, the 
magnitude of the high-momentum tails of both the mo- 
mentum distribution and the structure factor. In a trans- 
lation invariant situation, the contact is equal to n^ times 
the local value of the pair distribution function 5(0). Be- 
fore proceeding to derive these results from the OPE in 
explicit form, we note that our derivation remains un- 
changed if the sign of the interaction is reversed. All the 
results of this article can thus be extended to the case 
of an attractive Coulomb interaction by simply changing 



III. OPE FOR THE MANY-BODY 
WAVEFUNCTION 

The crucial physical insight, already implicit in the 
work of KimbalpSli^, relies on the intuition that the 
many-body wavefunction factorizes into a two-body con- 
tribution and a remainder whenever two particle coor- 
dinates come closer than the average interparticle dis- 
tance. In this limit, the two particles only feel their mu- 
tual Coulomb repulsion at short distance, with negligible 
corrections from the medium. This type of argument has 
in fact been used by various authora^SESl jj^ ^^j^g derivation 
of the Tan relations for Fermi gases with short range in- 
teractions. In order to prove the validity of this intuitive 
picture, we use the operator product expansion for the 
special case of the operator 



a(x)a(y)=V^(x)^i(y). 



(14) 



Inserting the general form (11) of the OPE, the A^- 



particle wavefunction ^ corresponding to an arbitrary 
TV-particle state \'^n) can be expanded as a formal power 
series 



VWNi 



1=3 



1 ^ 
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1=3 



where Wn (r) are the Wilson coefficients in an OPE of the 
operator '0^(x)V'^(y), which can be written as the sum 
of a spin singlet operator 

V'(x)v(y) = ^(VtW^;(y) - ^iWv^T(y)) (16) 

and a triplet operator which is symmetric in the spin in- 
dices. The leading order term in the OPE is associated 
with the operator ■0-0(0). Its Wilson coefficient can be 
obtained by taking the expectation value of Eq. ([16]) be- 
tween the vacuum and a two-particle state with (on-shell) 
energy p^ /m. The corresponding diagrams are sketched 
in Fig. [Tl We can express this matrix element as a one- 
body scattering wave function 



-*o 
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■*o 
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(a) 



(b) 



FIG. 1. (a) Contribution to the two-particle matrix element 
of the operator 0(— 1)7/)(|), denoted by the white circles, (b) 
Same for the operator ipip(0). The T-matrix is denoted by a 
gray rectangle. 



a result, it turns out that W^^(r) solves the two-particle 
s-wave Schrodinger equation 



-V^ 



1 
aor 



W^^{r) = 



(18) 



with boundary condition W^^(0) = 1. The solution of 



Eq. ( 18 ) up to linear order in r is 



W^^{r) 




(3D) 
(2D) 



(19) 



The OPE therefore provides a rigorous derivation of the 
intuitive short-distance factorization 



iTo*(-|'^=i'^; 



M/^^(r)*(0,t;0,;; 



(20) 



of the many-particle wavefunction ( [15| ) if two particles 
of opposite spin are close to each other. This factoriza- 
tion was previously considered by Lepage for two-particle 
systems in the context of effective field theorie^^ZJ, pj^e- 
cently, this was used to derive high-momentum tails for 
the unitary Fermi gas and the Coulomb gaa^S. 

Focusing on pure states, e.g., the ground state of jel- 
lium, the result ( 20 1 together with the fact that the two- 
particle wave function VF^^(r) is non-analytic at short 
distances implies power-law tails in the ground state mo- 
mentum distribution and the static structure factor. In- 
deed, for a pure state, the dimensionless and intensive 
momentum distribution is given by 



(0|V'(-|)V'(|)|p,t;-P,;) 

= (r|l + Gor|p,t;-p,;) = 



'0p(r), 



(17) 



using the bare retarded two-particle propagator Gq and 
the T-matrix T. In the second line, we have used the 
Lippmann-Schwinger equation for the scattering wave- 
function '0p(r) with energy p^ /m. Similarly, the matrix 
element of the operator tljip{0) is 0^(0)^^. We determine 
its Wilson coefficient by matching this matrix element to 
leading order in the energy p^ /m of the external state. As 



i(q) = ^ j d^R fd'^rfld'^n 



-zqr 



xvI/*(R_|,r2,...,r^)vI'(R-)-|,r2,...,rjv). (21) 

Its asymptotic behavior for large momentum is deter- 
mined by the integration regions in which both R + r/2 
and R — r/2 approach one of the particle coordinates r; 
for I = 2, . . . N simultaneously. By substituting the re- 
sult (201 in Eq. (21), the resulting high-momentum tail 



of the momentum distribution turns out to be given by 

n(q) 



ao 



q- 



'27ry n^g{0) 



ao 



q" 



(3D) 
(2D) 



(22) 



in accordance with the results in Refsi^^Hi^^'^SI 

The pair distribution function in first-quantized form 
reads (specifying to the translation invariant case) 



N 

g{r)^N{N-l) l[dW^{0,r,r3, . . . ^r^) 

•^ 1=3 



(23) 



Inserting Eq. 
function, Eq. 



9ir) = 



( |20[) in the definition of the pair correlation 
(1231), we obtain 







in r 
ao 27r 



5(0) 



(3D) 
(2D) 



(24) 



The pair distribution function at short distance there- 
fore exhibits a dip, rising linearly with slope g(0)/ao. 
For an attractive Coulomb force, where oq — ^ —a-Oi the 
dip is replaced by a cusp. The result (l24| agrees with 
the one obtained previously by Kimball''^ ''^^. Following 
Rajagopal et al^, we call this the cusp condition. It 
is interesting to note that a result which is completely 
analogous to that in Eq. (24) holds for Fermi systems 



with zero range interactions in one dimension, with the 
ID scattering length ai replacing the Bohr radius . 

The short-distance singularity in the pair distribution 
function leads to an asymptotic power law in the static 
structure factor 



S(q) = l + n f d'^re-''^" {g{r) - 1) 



(25) 



Substituting Eq. ( 24 ) , the static structure factor behaves 
like 



S{q) - 1 = 



Stt 715(0) 
ao g** 
in ng{0) 
ao q^ 



(3D) 

(2D) 



(26) 



at large momentum, where we have used that the Fourier 
transform of r in three and two dimensions is —Sn/q"^ and 
—2n/q'^, respectively, in the sense of distributions. The 
results in Eq. ( 26 1 are again in accordance with those 
obtained in RefiPhimai. 



The tails ( 22 ) and ( 26 ) in both the momentum distri- 



bution and the static structure factor are present for any 
state of the Coulomb system, a property that will be de- 
rived in detail in the following section. Here, we focus 
on the particular case of zero temperature and moderate 



values of r^., where the jellium ground state is a uni- 
form Fermi liquid. The characteristic momentum scale 
beyond which the asymptotic behavior applies is then 
set by the Fermi wave vector kp, which is related to the 
average interparticle distance rg via kp = l/(aro), with 
a = (4/97r)V3 ~ 0.521 in 3D and a = 1/V^ ~ 0.707 
in 2D. Introducing a dimcnsionless strength s of the tail 
in the momentum distribution via na-{q) — > s (kp/q)'^'^^'^ 
one obtainPil 



, , ■ „«8 5(0)r2 (3D) 
s(rs) = <( ^ 

2a^gi0y, (2D) 



(27) 



The dimensionless strength s is a continuous function 
of Ts, vanishing in both limits r<, — >■ and r<, — >■ oo. In- 
deed, the power-law tail in the momentum distribution 
is present even in the Wigner crystal, as long as the spin- 
polarization C remains less than one (recall Eq. (10)). 



The contact in this limit is expected to to vanish in an 
exponential manner with rs . In the opposite limit of 
high density, the fact that g{0) = 1/2 + 0(r,)Ei yields 
s ~ r^. The function s{rs), therefore, must have a maxi- 
mum, whose value appears to be much smaller than one. 
Indeed, according to recent quantum Monte Carlo calcu- 
lations of the momentum distribution in the Fermi liq- 
uid phase of the jellium model in 2E'^^' and 3E33!, ^j^g 
resulting dimensionless strengths s(10) ~ 0.006 (2D) or 
~ 0.009 (3D) of the power law tail are surprisingly small 
even at r^ = 10. In particular, they are almost two or- 
ders of magnitude smaller than the corresponding value 
s(oo) = 321n2/(37r2) ~ 0.749 of a Fermi gas with infinite 
short range repulsion in ID^^. 

In this context, it is important to mention that various 
approximate schemes which have been developed to treat 
the jellium problem at least in its Fermi liquid phase 
in fact violate the cusp condition (24). This is true in 



particular for the standard RPA approximation, which is 
exact in the long wavelength limit but violates the cusp 
condition, see e.g. Ref.l^. It is obeyed within extensions 
of the RPA which include local field corrections like the 
one developed by Singwi, Tosi, Land and Sjolandeil^^, 
at the expense, however, of violating the com press ibility 
sum rule S{q) ~ q^ at long wavelengths q — >■ CpE3. 

In contrast to zero range interactions, where the value 
of the contact determines the complete thermodynamics 
by a simple coupling constant integratiorP^^ the ground 
state energy of the jellium model requires knowledge of 
the pair distribution function at all distances. Interpo- 
lation schemes for the static structure factor and thus 
the complete pair distribution function g(r) which prop- 
erly account for both the long- and short-distance behav- 
ior of the homogeneous, unpolarized electron gas have 
been proposed by Gori-Giorgi, Sacchetti and BacheletP^ 
and may be used to develop improved versions of the 
exchange and c orrela tion energy functionals in density 
functional theorjl^^liSl 



IV. DIRECT OPE OF THE CORRELATORS 

In the following, we will show how the OPE can be 
used to perform an expansion of the one-particle den- 
sity matrix and the pair distribution function nh and ([8]) 
at the operator level. Apart from providing an alterna- 
tive derivation of the high-momentum tails which avoids 
discussing the many-particle wavefunction, this method 
makes evident a point stressed already in our introduc- 
tion: the short-distance properties derived here are valid 
completely independent of the state of the system. In 
particular, they hold in arbitrary few-body states or in 
equilibrium at any temperature. 

We start by considering the static structure factor, 
which - for q 7^ - is just the Fourier transform 



5(q) 



1 

N 



d^'R 



d^re-''5""( 



(R. 



R 



(28) 
of the density correlator. Equations ((sl) and (28) imply 
that its asymptotic behavior for large momentum q is 
dominated by the short distance behavior of the pair dis- 
tribution function. The pair distribution function is con- 
nected to the ti density correlator VJV'^ (~l) '^i^i (§) 
via definition (l8|. As shown in detail in App. BJ the 
short-distance OPE of this correlator to linear order in r 
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VIVKV-^W + . . . (3D) 



V'|VKV;(0) + . . . (2D) 



(29) 



where we have omitted the analytic term of order r since 
it does not contribute to the high-momentum asymp- 
totics. The OPE ( 29 ) , together with the definitions of the 
pair correlation function Q and the static structure fac- 
tor ( 28 ) , reproduces the high-momentum behavior ( 26 1 . 



In particular, when taking the expectation value of (29 1, 
the contact operator i/j^iplip^ipJO) produces the pair cor- 
relation function at zero separation. 

The momentum distribution n^- (q) describes the prob- 
ability to find a particle of spin a with momentum q. In 
second quantization, it is defined as the Fourier transform 
of the one-particle density matrix: 

n.(q) ^i I df^x j d^ye-^^-^y--^ {l^'H^^v)) 



V 



^Jd^Rjd''re-^^'^{i;i (r - |) ^, (r - 



(30) 



The non-analytic Wilson coefficients in a short-distance 
expansion of the one-particle density matrix therefore de- 
termine the high momentum tail of n^ (q) . A quite elab- 
orate calculation, which is discussed in detail in App. [B| 
shows that the OPE of the one-particle density matrix is 
given by 



*i(-D*'(D 



e-§-^Vj(0)] [e5-^V,(0) 
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(3D) 



(31) 



V'jV'jV'tV';(0) + . . . (2D) 



The first Wilson coefficients of the bilinear operators are 
the coefficients in a Taylor expansion of the operators on 
the left-hand side. The contact operator Oc defined in 
Eq. ( 13 ) is the leading order term associated with a non- 
analytic Wilson coefficient of order 0{r^) and 0{r'^ logr) 
in 3D and 2D, respectively. Substituting Eq. (31 1 in (30 1, 
we obtain 



The fact that the non-analytic terms in Eq. (31) appear 



na{(\) = 



'47r\ n^n;ff^;(0) 
an) q^ 

'27ry n^n;5^;(0) 



-,6 



(3D) 



(2D) 



(32) 



for an arbitrary state with a possible non- vanishing spin 
polarization ^. Summing over a —'[, J,, we recover the 
previous result for the momentum distribution at large 
momentum in the spin-balanced Coulomb gas, Eq. (22). 



at the level of operators shows that the tails in the mo- 
mentum distribution are also present in phases where 
translation invariance is broken, in a few-body situation, 
or at arbitrary temperatures. In a Wigner crystal, for 
instance, the product n^n^g-\i{0) has to be replaced by 

iy'd^i?(V'|^|^^^;(R)), (33) 

which is again an intensive quantity in the thermody- 
namic limit N,V ^ oo at fixed average densities n^,n^. 
In a few-body situation, in turn, these densities vanish 
but there is still a finite expectation value of the contact 
operator. For the hydrogen atom for instance, one finds 



i?(V'jV'|V'tV'i(R)) = l^n00( 



(34) 



in agreement with the result derived in the introduction. 



V. HIGH-TEMPERATURE VERSUS 
CLASSICAL LIMIT 

Beyond the derivation of exact relations which con- 
strain the short-distance properties of Coulomb systems 
in quite general terms and which - as has been shown in 
the preceding sections - all involve the contact (OcCR.)), 
quantitative results for specific phases of jellium or non- 
trivial few-body Coulomb systems require to calculate 
the value of the contact as a function of both the in- 
teraction strength rg and temperature T. Since relevant 
values of r^ are beyond the regime where perturbation 
theory can be applied, this can only be achieved numeri- 
cally, for instance via quantum Monte Carlo calculations, 
see e.g. Refs.l^SlMl fgj. gome recent results. In the follow- 
ing, we calculate the value of the contact in the classical 
and the high-temperature limit. Surprisingly, it turns 
out that for Coulomb interactions these two limits are 
not equivalent. In fact, they turn out to be completely 
opposite. 

Consider the Coulomb gas in the regime 



/cbT> 



h^n^/^ 



(35) 



where the thermal energy is much larger than the de- 
generacy energy. This is the standard limit of a non- 
degenerate gas, in which the average interparticle spac- 
ing n"^/'* is much larger than the thermal wavelength 
At = ;i(27r/mfcsr)i/2: 



i^/'^Xt < 1. 



(36) 



In this limit, thermodynamic properties can be calculated 
by expanding in powers of the fugacity z = exp (/3/i) = 
nX^/2 ^ 133. The non-degeneracy condition (36) does 
not involve the strength e^ of the interaction and is satis- 
fied both in the infinite temperature and in the classical 
limit. Now, for systems with Coulomb interactions, there 
is a second and purely classical, so-called Bjerrum length 
£b = e^ /{ksT), which - keeping h finite - eventually be- 
comes shorter than the thermal length at sufficiently high 
temperatures. As a result - already noted in Refl^Sl- the 
order in which the limits T — )> c» or S — >■ is taken mat- 
ters. Taking T — ?> oo before /i — >■ results in the following 
hierarchy of length scales: 



I ^/"^ > ao > At > ^i 



(37) 



In turn, taking the classical limit h -^ before T ^- oo, 
we find: 



-l/d 



> ^B > At > ao- 



(38) 



As will be shown below, these two limits give quite dif- 
ferent results for the value of the contact. Since the gas 
is non-degenerate in both cases, the contact value of the 



pair distribution function can be calculated to leading 
order in the virial expansion, which just involves an in- 
tegration 



n^g{0) 



_ 2^ f^^ ,^./„ 



\rp 



{2TThy 



iVp(0)|2, (39) 



o;3m 



of the square of the relative Coulomb wavefunction 
|V'p(0)P at the origin with the classical Boltzmann dis- 
tribution for the relative momentum p. Here, z 
the fugacity while |?/'p(0)p is given by^ 

2tt^ 

l^p(0)P 



e27r5 _ I 

2 
, e^-^i + 1 



(3D) 

f2D) 



(40) 



with ^ = me^ /2hp. In the infinite temperature limit, the 
integration gives 



ff(0) 




V27T- 



1- 



V2h 



(3D) 
(2D) 



(41) 



which is the expected result for a classical ideal gas. The 
corrections involve the ratio Ib/Xt and thus vanish for 
large temperatures like ^ 1/\/T. Note that, although 
g{0) is of order one in this limit, the high-momentum 
tails in this limit are present only for very large q ^ X^ . 
More relevant for low temperature, non-degenerate 
plasmas is the classical limit fi. — >■ for which Ry^ fesT. 
This limit is reached, for example, for Coulomb gases of 
charged dust particles in astrophysics^. In this limit the 



integration in (39) gives rise to a contact which vanishes 



exponentially like 



m = <^ 



31/2 l^A^ 
27r2i/3 f£B_ 
At 



31/2 



4/3 



1/3 



3-R 
'21/3 



37r 
-31/3 



2/3 



2/3 



(3D) 

(2D) 
(42) 



Equations (41 ) and (42) are the main results of this sec- 
tion. The Wilson coefficients of the momentum distri- 
bution and the static structure factor diverge as l/a§ ^ 
1/h'^ and l/flQ ~ 1/h^ as fi — >■ 0. The high momentum 
tails are expected to occur for q ^ X^ . Thus, their char- 
acteristic momentum scale is pushed to infinity as h —> 0. 
The total weight of the tails, containing both the Wil- 
son coefficients and 5(0), is exponentially suppressed and 
guarantees a well defined kinetic energy in this limit. We 
point out that different high-momentum tails are present 
in the classical system. Indeed, the structure factor of 
the classical electron gas decreases as 1/g^ at intermedi- 
ate momentum k < q < 1/Ib, where k = \/nlB is the 
inverse Debye-Hiickel lengtHSl. 



A. Diagrammatic derivation 

In the infinite temperature limit, the virial expansion 
of the contact g{0) in Eq. (41 1 can also be obtained using 
a diagrammatic formalism. This method was originally 
introduced by Vedenov and Larkin to derive the equa- 
tion of state of an electron gas^, and was recently used 
by various groups to determine the virial exp ansion of a 
Fermi gas with short-range interactionfpfiHm^ 

The perturbation series of an arbitrary correlator in- 
volves all Feynman diagrams that connect to the operator 
insertions. In the absence of a small parameter, this gives 
rise to a very large number of diagrams. The key point of 
the method is that in the infinite temperature limit, this 
number is drastically reduced by exploiting the causal 
structure of the propagators, which are defined as 



G(T,k) = 



-(1 - nk)e-(^''-^)^ 



n^e 



-(ek-M) 



T> 
T< 



(43) 



where rik = l/{e^^'^^~'^^ + 1) is the Fermi-Dirac distribu- 
tion. In imaginary time, they carry the whole dependence 
on the fugacity and, thus, an expansion of a diagram in 
the fugacity corresponds to an expansion of the propaga- 
tors: 



G{t, k) = G(°) (r, k) + G<-'Ht, k) + 0(z2), (44) 



with 



G(°)(^,k) = -e(T), 
GW(r,k) = 



and 



ze''^e-''''^'^+^\ 



(45) 
(46) 



Q is the Heaviside step function. The coefficients G^^' 
with i = 1, 2, . . . can be treated as separate diagrammatic 
elements, e.g. by denoting them by a line that is slashed 
i times. The order of a diagram in z is then set by the 
sum J2i ^i I where iV^ is the number of propagators of 
type i. It is important to note that the leading order 



term ( 45 1 is purely retarded, and while the calculation 



of a general order may be unwieldy, only a very limited 
number of diagrams contribute to the virial expansion of 
a correlator to leading order in z. 

Figure [2] shows the two leading order contribution in 
e^ to the pair correlation function 



m(0) = ^- {Trn^{P)n^m 



(47) 



where we define the density operator in the usual sense 
7i(j = — lim^^o- Tripa {T)ij}t{0) to avoid an ordering am- 
biguity—. Figure 2(a) corresponds to the noninteracting 
result 

5n(0) = 1' (48) 

where we use the relation 71^ — z^yX^ . The 0{e^) con- 





FIG. 2. (a) 0{e°) and (b) O(e^) contribution to g{0). Imag- 
inary time runs from the left to the right. Black dots denote 
the density operator, dashed lines the instantaneous Coulomb 
interaction, and (slashed) lines the particle propagator as ex- 
plained in the text. 



tribution in Fig. |2(b)] reads 

aim - 



V2n^ (3D) 



V2Xt 



(49) 



(2D) 



which coincides with our previous result (41). Note that 



higher order contributions to g{0) contain infrared diver- 
gences. They can be removed by summing the divergent 
parts of an infinite number of ring diagrams^, which 
gives rise to Debye-Hiickel corrections that are of higher 
order in the density. 



VI. SUMMARY AND CONCLUSIONS 

In summary, we used a short-distance operator prod- 
uct expansion to derive the high-momentum tail of the 
structure factor and the momentum distribution of the 
Coulomb gas. Since these results are based on operator 
identities, they hold for pure states as well as for mix- 
tures, and, in particular, for different phases, such as a 
Fermi liquid or a Wigner crystal. The key idea behind the 
derivation of our exact results is the separation of short 
and long distance scales. The functional dependence on 
the high-energy scales can be calculated exactly while 
the low-energy contribution factorizes. This multiplica- 
tive constant is the contact value of the pair distribution 
function, consistent with the intuitive expectation that 
the short-distance physics is determined by the probabil- 
ity to find two particles at the same point. These results 
are in close analogy to the Tan relations for zero range 
interaction&H', in which an analogous contact enters the 
coefficients of the high-momentum tails. 

Furthermore, we calculated the contact in a virial ex- 
pansion. It turns out that there are two possible limits 
which yield different results: on the one hand, in the 
high-temperature limit, the contact approaches the ideal 
Fermi gas value with power-law corrections in the tem- 
perature. In the classical limit, on the other hand, the 



contact vanishes exponentially, which is crucial to en- 
sure a well-defined transition to the classical regime in 
which the Coulomb repulsion between the particles pre- 
vents them from being at coincident points, and thus no 
short-range correlations are present. 

The universal relations obtained in the present article 
are by far not exhaustive. Indeed, many more relations 
could be derived within the framework introduced here. 
A short-time OPE analogous to Refs.l^SHSl] -^^ould give re- 
sults for dynamical correlators, such as, for example, the 
current response function or the dynamic structure fac- 
tor, which display short-range correlations that are not 
captured in a random phase approximatioipSl. A simi- 
lar analysis can be carried out for the spectral fun ction, 
which possesses a high-frequency tail as derived irP^ESl. 
Beyond applying the OPE to the Coulomb gas, it should 
be straightforward to generalize the results in this article 
to other many-fermion systems with long-range interac- 
tions, such as quantum gases of dipolar particles, which 
can be studied experimentally^SHlIl 
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Appendix A: Feynman Rules 

This appendix summarizes the Feynman rules of the 
Hamiltonian m in momentum space. An energy w and 
a momentum q are assigned to each internal and external 
line. The bare propagator is denoted by a continuous line 
and contributes a factor Go(w, q) = i/{uj — q^ /2m + iO). 
We represent the interaction between two fermions by 
a dashed line. It contributes a factor Airie'^/q'^ in 3D 
and 27rie^/g in 2D, where q is the difference between 
the center of mass momenta of the ingoing and outgoing 
fermions. Finally, each undetermined momentum and 
energy is integrated with measure / d'^q/ {2ttY J duj/2TT. 

The T-matrix insertion iT{p,p' ,k) is denoted by a 
gray rectangle, where p and p' are the center of mass 
momentum of the ingoing and outgoing atoms, respec- 
tively. We denote the center of mass energy of the ingo- 
ing atoms hy E = k'^/m. Note that this energy is not 



necessarily on-shell, i.e. we do not impose the condition 
p'^. The T-matrix solves the Bethe-Salpeter 



p^ 



equation, depicted diagrammatically in Fig. 3(a) and 

j.gg^(j^25|26l59 60 



zT(p,p',fc) ^ 

Jo dxH{p,p',k) 






-Aiire^k / dxx"^^ 



X 



1/2 



dx m/'^{p.,p\k) 



(3D) 



(2D) 



(Al) 



where i?(p,p',fc) :^ ^Pip+p'fx-ik"^ -p^){k^ -p'^){l- 
x)^ and f = vie^ /2k is called the Sommerfeld parameter. 



Appendix B: Operator Product Expansion 

In this appendix, we collect the matrix elements and 
momentum integrals needed to perform the OPE for the 
Coulomb gas. For simplicity, we set fi = 1 in the fol- 
lowing. Since the operator product expansion is state- 
independent, it is sufficient to evaluate the matrix ele- 
ments of the operator products ^ and (IS]) and the local 
operators between selected few-particle states. The Wil- 
son coefficients are determined by matching the terms in 
an expansion of these expectation values in the external 
parameters of the state. The Wilson coefficients of the 
bilinear operators, i.e., operators that contain one field 
operator and its hermitian conjugate, are the coefficients 
in a Taylor expansion of the operator product, which can 
be obtained by matching the matrix elements between 
a one-particle state. To compute the contact's Wilson 
coefficient, we choose a two-particle state with zero rel- 
ative momentum and (off-shell) energy k'^ /ra^ which we 
denote by {k\0\k). These matrix elements are the sum 
of four diagrams with either scattering or no scattering 
in the initial and final states. Some of them are depicted 
in Fig. [3} 



a. Structure factor 



As explained in Sec. IV the short-distance behavior 



of the pair correlation function is inferred from an OPE 
of the ti density correlator n^(— |)n^(|). The matrix 
element of this operator between a two-particle state can 
be expressed in terms of diagram |3(b)[ The complete 
matrix element is the sum of four diagrams with either 
scattering or no scattering in the initial and final states: 



{k\n^{- 



■l>H{\m 



duj 

2^ 



d'^q 

(27r)° 



e'^'^iTiO, q, k)Goiuj, q)Go{E - u, -q) 
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^o*- 



^o*- 





-s^ 





*o*~ 



(a) 



(b) 



(c) 



(d) 



(e) 



FIG. 3. (a) Bethe-Salpeter equation for the T-matrix, which is indicated by a gray rectangle, (b-e) Contribution to the two- 
particle matrix elements of (b) the density-density correlator nf{—^)n_i{^), (c) the contact of the pair distribution function 
g{0), (d) the one-particle density matrix ipt{—^)i'cr{^), and (e) bilinear operators, respectively. 




a-oj 


[4/ dxx^^'^ ,^ "^ 
Jo dx{l + x 


)\ 


+ 0(r2) 


(3D) 


2r\ 
aoj 


Jo dx\ + x^ 


+ 0{^) 


(2D) 



(Bl) 



where we used the integrals 



and 



i 1 



d^q e^i"" 

(27r)3 (g2-a2)(g2_52) " 4^ a + 6 " 8^ 



+ 0(r2) 



(B2) 



(27r)2 (g2_a2)l/2(^2_52) - ^ ^^2 _ ^2 " 2^ 



+ 0(r2). 



(B3) 



The factor in square brackets in Eq. (Bl ) depends on the details of the states and must not contribute to the Wilson 



coefficients. It is matched by the expectation value of the contact tp^.tplip^iliAQ). Consider the diagram in Fig. 3(c' 



1 x< 



— 4i^ / dx , , „ 

{l + xf 



J (27r)'' J 27r I o-^ /' j 2; 

-2i^ / dx 



l + x 



(3D) 

(2D) 



(B4) 



The full matrix element of f/'If/'IV'tV'i (0) contains three additional diagrams: 



(fc|V'|v|V'tV';(0)|fc) 



where we have expanded the result to leading order in k. 



Momentum distribution 



hf'dxx^^ '^ ,^ "^ 
Jo dx (l + x 


r. 


^ 1 


(3D) 


I f^ : d a;i/2 1 
2 / dxx'i 7 
Jo dxl-l-a; 


2 


1 


(2D) 



(B5) 



To obtain the asymptotic form of the momentum distribution (30), one performs an OPE of the nonlocal operator 
'0j-(— f )V'(T(f ), whose expectation value gives the one-particle density matrix 7(|,— |), cf. Eq. (It]). Since insertions 
of this operator on external legs are matched by bilinear operators, the only relevant diagram that contributes to the 
Wilson coefficient of the contact operator involves scattering in both initial and final state as shown in Fig. |3(d)[ As 
we are only interested in the leading order non-analyticity of the Wilson coefficient of the zero distance pair correlator, 
we expand the T-matrix, as well as our diagram, as a power series in k around A: = 0. This procedure introduces 
infrared divergences, which we regulate by introducing an infrared cutoff /u. The expansion of the T-matrix with 
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X 0-» 1 K X 1 — »— O K 

I I 

X 0-» — I X X 1 — »— o K >- 

(a) (b) 



X X 



(c) 



FIG. 4. Leading-order contributions to the structure factor (a,b) and momentum distribution (c) 



respect to fc ~ 1/^ is given by 



27re2 ,. 

I ^„ ^ — lUTL t- 



27re^ 



lim i- 



ne 



ne 



2eq M^O ^2^g2+^2 



(3D) 
(2D). 



(B6) 



In addition, we expand our propagators as 1/(9^ — fc^) = l/((7^ + /i^) + 0{k^). In the limit k 
given by 



0, diagram 3(d) 



IS 



(2^ 



duj 
2n 



e^i-"-(*r(0, q, fc ^ 0)fGoiio, q)^GoiE - uj, -q) 



480^^4 



.51 



leage 



r r 
8^1 " 15 
1 r2 r4 /3 



2^3 



2/i4 8Ai2 32 V 4 



0(r*^ 
-log 



(3D) 



e'^'^ fir 



(B7) 



-0(r5) (2D). 



The matrix elements of the one-particle operators (Fig. |3(e) I match the analytic terms in this expansion. The 
remainder is of order 0{r^) and ©(r^logr) and is matched by the contact operator: 



WoSr) = 



47ry r^ 



ao J 28807r 
2ttV r* 



+ 0(7-6) 



Oo 



1287r 



log- 



pjE It 



(3D) 



-0{r^) (2D) 



(B8) 



c. Green's function OPE 

In addition to the derivation outlined in the previous 
section, the high-momentum tail of the structure factor 
can also be obtained from a short-time and -distance 
OPE of the time-ordered density Green's function 



iG„(w,q) = I dt I d'^3 



gicjt — -iq-x 



(rn(t,x)n(0,O)). 



(B9) 



For q 7^ 0, it is related to the structure factor by 



5(q) - 1 = - hm 
n t-i-o- 



diu 

2n' 



zG„(w,q). (BIO) 



The integral is evaluated by closing the contour in a large 
semicircle in the lower half of the complex w-plane. Only 
Wilson coefficients with poles in both half-planes con- 



tribute to the high-momentum tail. In the limit of negli- 
gible external scales, diagrams have vanishing residue if 
they can be traversed from one operator insertion to the 
other by following the fermion lines. The Wilson coef- 
ficient of the operator Oc are read off directly from the 
diagrams in Figs. 4(a) and|4(b)] 



WoA^i^. 




1 

iO){-uj- 
1 



iO) 



q^ (cj - Eq + iO)(-w - Eq + iO) 



(3D) 

(2D) 



(Bll) 



where e^ = q^ /2m, and the external lines couple to the 
contact (7-|-|(0). Performing the contour integration re- 
produces the result ( 26 ) . 



We can apply a similar argument to determine the 
high-momentum tail of the momentum distribution. 
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which is related to the single-particle Green's function 

*G,(i,x) = (rV..(i,x)V;t(0,O))bjEl 



n,(q)--lini / ^ e^^^* *G,(q, w). (B12) 



This relation was used to derive the momentum distribu- 
tion of a Fermi gas with short-range interaction^^. The 
first nonzero contribution is given by the contact opera- 



tor, which has the Wilson coefhcient (cf. Fig. 4(c) ) 
W^o<:(w,q) 



-1 



47rp/ 

g2 ) (tu-£q-HiO)2(-a;-eq-l-iO) 
'27rp/ 



(w - Eq + i0)2(-a; - £q + iO) 



(3D) 

(2D) 



(B13) 



Calculating the residue in Eq. (B12) yields the previous 
result (|32l). 
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